clear; clc; close all;

% 参数设置
A_values = [0.05, 0.1, 0.2];  % 电流源幅值A
omega_values = [0.1, 0.2, 0.4];  % 电流源频率omega
n = 1:500;  % 迭代步数
q0 = 0;  % 初始状态
iterations = length(n);  % 迭代次数

% 初始化内在状态矩阵
q = zeros(length(A_values), iterations);  % 存储内在状态
v = zeros(length(A_values), iterations);  % 存储输出电压

% 创建图形
figure;

% 子图(a) 电流源幅值变化
subplot(1, 2, 1);
hold on;

for i = 1:length(A_values)
    A = A_values(i);
    q(i, 1) = q0;  % 初始状态

    % 计算内在状态q和输出电压v
    for j = 2:iterations
        i_n = A * sin(omega_values(2) * n(j));  % 电流源i_n（固定频率omega=0.2）
        q(i, j) = q(i, j-1) - i_n;  % 更新内在状态
        v(i, j) = cos(q(i, j)) * i_n;  % 计算输出电压v
    end
    
    % 绘制存储器特性曲线
    plot(q(i, :), v(i, :), 'LineWidth', 2);  % 绘制每个A值的曲线
end

% 设置图表标题和标签
title('Memristor Hysteresis for Different A');
xlabel('q_n');
ylabel('v_n');
legend(arrayfun(@(A) sprintf('A = %.2f', A), A_values, 'UniformOutput', false), 'Location', 'best');
grid on;
hold off;

% 子图(b) 电流源频率变化
subplot(1, 2, 2);
hold on;

for i = 1:length(omega_values)
    omega = omega_values(i);
    q(i, 1) = q0;  % 重置q

    % 计算内在状态q和输出电压v
    for j = 2:iterations
        i_n = A_values(2) * sin(omega * n(j));  % 使用A=0.1作为示例
        q(i, j) = q(i, j-1) - i_n;  % 更新内在状态
        v(i, j) = cos(q(i, j)) * i_n;  % 计算输出电压v
    end
    
    % 绘制存储器特性曲线
    plot(q(i, :), v(i, :), 'LineWidth', 2);  % 绘制每个omega值的曲线
end

% 设置图表标题和标签
title('Memristor Hysteresis for Different \omega');
xlabel('q_n');
ylabel('v_n');
legend(arrayfun(@(omega) sprintf('\\omega = %.2f', omega), omega_values, 'UniformOutput', false), 'Location', 'best');
grid on;
hold off;
